%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PNJunction() is a function that computes the charge and potential
% distribution of a p-n junction formed by two arbitrary materials
% Inputs : 
%   Q1, Q2 : Total background charge density of material (nuclei + dopant)
%   DOS1, DOS2 : DOS defined on energy grid for each material
%   KT : Temperature
%   mf : potential mixing factor
%   N : Spacial Grid Size
%   L : Periodic Length
% Outputs :
%   V : Potential Distribution
%   R : Charge Distribution
%   Err : Record of error in each iteration

function [V,R,Err] = PNJunction(Q1, Q2, DOS1, DOS2, KT, mf, N, L, varargin)
%%%%%%%%%% Function Setup %%%%%%%%%%
%% Geometric Setup (Real Space Grid, Fourier Grid)
if mod(N,2) == 1             % Making the size N an even number to make the k 
    N = N - 1;               % grid symmetric about the origin
end
% Fourier Grid
% Note here that for an even Fourier Grid, the DC frequency component is at
% the position 1+N/2.
n = -(N/2):(N/2-1);
oddk_index = 2:2:N;        % index of the odd k points
evenk_index = 1:2:N-1;         % index of the even k points
kn = n*(2*pi/L);
% Real Space Grid
h = L/N;                     % grid size in real space (in Bohr Radius)
x = linspace(-L/2, L/2,N);   % Real Space Grid
nc = 1+N/2;                    % The middle index 
%% Energy Grid Setup/DOS Setup
if ~isequal(size(DOS1), size(DOS2))
    error('Energy Grid:', 'The two DOS have different size!');
end
E0 = DOS1(1,:);
D1 = DOS1(2,:);
D2 = DOS2(2,:);
E = linspace(E0(1),E0(end),1E5);
% Pre-computing the surface charge as a function of an ultra-fine energy
% grid, stored as S1,S2
tic;
S1 = BandCharge(D1,E0,E,KT);
toc;
S2 = BandCharge(D2,E0,E,KT);


BGCharge = zeros(1,N);       % The total background charge
BGCharge(1:nc-1) = Q1;
BGCharge(nc:end) = Q2;
Q1BG = Q1*ones(1,N); 
Q2BG = Q2*ones(1,N); 
Ef1 = FindFermi(x,zeros(1,N),Q1BG,E,S1,S1,0);
Ef2 = FindFermi(x,zeros(1,N),Q2BG,E,S2,S2,0);
display(Ef1);
display(Ef2);

%% Physical Constant Setup (In Hartree Atomic Unit)
eps = 1;            % Relative permittivity
eps_0 = 1/(4*pi);   % Vacuum permittivity 
e = 1; % Electronic Charge
hbar = 1;
D2Kern = 1/(2*eps*eps_0);
%% Initial Setup for v0 in k space
%% Square Wave in k space (for symmetry)
if isempty(varargin)
    qn = atan(x)/(0.5*pi);
    qn(1:N/4) = atan(-(x(1:N/4)+L/2))/(0.5*pi);
    qn(3*N/4+1:end) = atan(-(x(3*N/4+1:end)-L/2))/(0.5*pi);
    v0n = qn*(Ef2 - Ef1)/2/e;
    v0 = fftshift(fft(fftshift(v0n)))/N;
    %figure();
    %stem(imag(q));
    %plot(N*real(fftshift(ifft(ifftshift(v0)))));
else
    vin = varargin{1};
    if ~isequal(length(vin), N)
        error('Spatial Grid :', 'Input voltage has wrong size');
    end
    vx0 = varargin{1};
    v0 = (1/N)*fftshift(fft(fftshift(vx0)));
    
end
Ef = 0;
% Error 
err = Inf;
Err=[];
%%%%%%%%%%%%%%%%%%%%%%%% End of Setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Function Body
% Setting up figure for debug
figure;
ax1 = axes('box','on','XTickMode', 'Manual','XTickLabel',{},'xlim',[1,N],'Position', [0.15,0.55,0.8,0.35],'Visible', 'on');
ax2 = axes('box','on','xlim',[1,N],'Position', [0.15,0.10,0.8,0.35],'Visible', 'on');

it = 0;
while err>1.0e-7
%for i = 1:300
    it = it + 1;
    vx = N*real(fftshift(ifft(ifftshift(v0))));
    % v-> DOS -> s
    % Determine the fermi level, obtain the adjusted charge distribution
    Ef = FindFermi(x,vx,BGCharge,E,S1,S2,Ef);
    sx = BGCharge-V2Charge(S1,S2,E,vx+Ef);

    % Fourier Transform of charge (on x)
    sk = (1/N)*fftshift(fft(fftshift(sx)));
    % s -> Gauss' Theorem -> v
    v = D2Kern*sk./abs(kn);
    v(nc) = 0;             % Setting the center frequency to 0 (no DC term)

    % Recording the error
    err = sum(abs(v-v0))*h;
    Err = [Err,err];
    
    % mixing
    v0 = (1-mf)*v0 + mf*v;
    
    % Display Result for Each Iteration
    %display(sprintf('Ef(%d)=%.15f err(%d)=%.15f maxV %e: maxv0: %e',it,Ef,it,err,max(V),max(v0)));
    if mod(it,100)==0
        display(sprintf('Ef(%d)=%.15f err(%d)=%.15f',it,Ef,it,err));
        MyDraw(vx,sx,N,Ef,Ef2,Ef1, ax1, ax2);
    end

end

% Outputs
V = N*ifftshift(ifft(ifftshift(v0)));
R = N*ifftshift(ifft(ifftshift(sk)));



end


%% My Functions
%%%%%%%%%% Function MyDraw %%%%%%%%%%
function MyDraw(Vx,Sx,N,Ef,Ef2,Ef1, Ax1, Ax2)
    Vmax = (Ef2-Ef1)/2; 
    nc = N/2;
    axes(Ax1);
    plot(Ax1, real(Vx)+Ef); hold(Ax1, 'on');
    plot(Ax1, ones(1,N)*Vmax+Ef,'--k');
    plot(Ax1, -ones(1,N)*Vmax+Ef,'--k');
    plot(Ax1, zeros(1,N), '--k'); 
    ylabel('V(x)'); hold(Ax1, 'off');
    %xlim([nc-100,nc+100]);
    
    axes(Ax2);
    plot(Ax2, real(Sx));
    ylabel('\sigma(x)');
    %xlim([nc-100,nc+100]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


